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I. 


INTRODUCTION 


The  Discrete  Convolution  Method  (DCM)  is  an  iterative 

solution  technique  for  solving  the  matrix  equation 

formulated  by  using  the  Method  of  Moments  (MOM)  [1],  This 

is  accomplished  essentially  by  looking  at  the  matrix 

equation,  as  a  (set  of)  convolution  equation(  s)  .  The  DCM 

can  solve  a  properly  formulated  NxN  matrix  equation,  with 

3 

only  NlogN  order  multiplicative  operations,  instead  of  N  as 
in  Gaussian  Elimination.  The  number  of  iterations  needed 
for  a  given  accuracy  is  also  found  to  be  practically 
independent  of  the  size  N. 


II.  FORMULATION 

We  first  prove  that  a  properly  formulated  MOM  matrix 
equation  can  be  viewed  as  a  convolution  process  and  then 
develop  the  solution  technique.  Depending  on  the  type  of 
problem,  the  equation  can  be  reformulated  into  two  types  of 
convolution  processes. 

1)  For  some  MOM  problems,  the  matrix  equation  can  be 
rewritten  as 
N 

E  Zmn  Jn  =  Vm  ,m=1  »2*  •  • • *N  (1 ) 

n=  1 

where  V  '  s  are  all  known.  For  most  of  these  problems  it  is 


v 


il'wT Ki VTV5  **.  *."  •-  . 

** • - .  *  » 
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1  possible  to  choose  expansion  functions,  renumber 

them 

and 

I  add  dummy  segments  (if  needed),  so  that  (1)  becomes 

1  "l  N2 

nm 

t.  n  •• 

1  V1  V1 

•  V  z  „  j 

n  =1  P1P2*  *  *PMqiq2‘ *  *qM  qiq2“*qH 

K 

Vf>1P2*  *  *PM 

(2) 

|  where  p1 =1 , 2 , . . 

.  ,H1 

|  ^2=^*^*** 

.  ,n2 

|  •  • 

I  p..=  1  ,  2  ,  .  . 

•  ,NU 

•  H 

|  Here  p1p2...pM  are  a  renumbering  of  the  original  N 

segments 

jg  plus  the  added 

dummy  segments,  a  total  of  N  ^ 

N2  *  *  • 

V 

|  Furthermore 

|  Zp1P2***PMq 

1q2’**qM  Pl"qi ,P2-q2’ ’ * *PM~qM 

(3) 

|  and  V 

are  not  all  known.  The  values 

of 

V' s 

I  P1 P2  *  *  *  PM 

1  corresponding 

to  dummy  elements  are  unknown 

but 

J '  s 

I  corresponding  to 

these  are  zero.  If  we  call  the  domain 

of 

I  original  segments  S  and  the  domain  of  all  segments  S 

■  ,  then 

I  v 

are  known 

P1 P2  *  ’  '  PM 

JplP2*  *  *PM 

are  unknown  (to  be  solved  for) 

(4) 

if  p 1 p2  *  * ’ PM C  S 

and 

V 

are  unknown 

P1P2*  *  *PH 

Jpl P2 ’ *  *  PM 

are  known  (sO) 
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if  p,p_. . .p„e  (s  -s) . 

l  c.  n  e 

Combining  (2)  and  (3).  we  get 
M1  *2  MM 

H  XI  -  H  z(prqi,p2“q2,,*,ptrqM)J(qi,q2,**-qM) 
n1 =1  n2=1  nM=l 

:  ^ ( p^ ,Pj i • • • i P^)  (5) 

where  p^  =  1  ^....N^pgsl  ,  2 , .  .  .H2;  ...pM=1  ,2  .... 

But  (5)  can  be  easily  recognized  as  an  H  dimensional 
convolution  equation  [2].  Therefore  (5)  can  be  rewritten  as 

1  «  7  =  “v  Plp2...pMcse  (6) 

where  "  denotes  convolution  of  the  appropriate  (M  here) 
order  or  dimension.  To  give  a  one  dimensional  example, 
consider  scattering  from  two  quarter  wavelength  straight 
wires  that  lie  along  a  single  axis  0.15  wavelength  apart,  as 
shown  in  Fig  .  1(a). 


(gap) 

Fig.  1(a)  Two  quarter  wavelength  wires  seperated 


by  0.15  wavelength  gap 


»  I  I  1 


I1  I  I 


%  T  * - 
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Fig.  1(b)  Segmentation  used  for  the  problem 
of  Fig.  1(a) 


Further,  suppose  that  we  break  the  wire  into  0.05 
wavelength  segments  and  also  add  three  dummy  segments  in 
between  as  shown  in  Fig.  1(b).  The  matrix  equation  is  then 


T  ZmnJn  *  V  m  .»=1  *2 . ^ 

£—4  mn  n  m 


However,  it  is  clear  that  the  value  of  Zmn  depends  only  on 
(m-n) .  Therefore  (7)  can  be  rewritten  as 


£  Z(m-n)  J  (  n)  =  V(m)  ,m=1,2 . 13 


In  other  words,  the  equation  is  a  one  dimensional 
oonvol ution , 

Z  •  J  s  V  (9) 

where 

V  are  known 

tn 

J_  are  unknown  for  m=1  , 2 , 3 .  ,5 ,9 , 10 , 1 1  , 12 , 13 
hi 

V_  are  unknown 
m 

Jn  are  known  for  m=6,7,8  (10) 

Other  examples  of  the  one  dimensional  convolution  are 
helical  wires,  infinite  strips,  infinite  circular 
cylindrical  segments,  linear  antenna  arrays,  etc. 

A  two  dimensional  example  is  the  MOM  formulation  of  a 
linear  antenna  array  problem  using  several  expansion 
functions  for  each  antenna,  as  shown  in  Fig.  2(b). 


Fig.  2(a)  A  linear  array  of  wire  dipoles 


ms* 


v'v'\\vv 
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Fig.  2(b) 


Expansion  functions  used  for  the  antenna  current 


If  p1  denotes 
function  number  for 


p1p2q1q2 


antenna  number 
each  antenna, 
Z< Pt -q i ,P2-q2> 


and  denotes  expansion 

it  is  easy  to  see  that 
(11) 


Note  that  if  we  use  only  one  expansion  function  per  segment 
we  get  a  one  dimensional  convolution  equation. 


A  three  dimensional  example  is  the  solution  of  a 
rectangular  antenna  array  problem  by  using  several  functions 
per  antenna.  The  problem  is  as  shown  in  Fig.  3. 


-•  -r 


U- 


Fig.  3  A  planar  array  of  wire  dipole  antennas 

Mote  that  if  we  use  a  single  expansion  function  per  antenna 
we  get  a  two  dimensional  convolution  equation.  A  four 
dimensional  equation  would  be  produced  by  cubic  array 
problems  using  multiple  expansions,  etc. 

2)  Most  of  the  remaining  MOM  problems  can  be  reformulated 


For  the  first  example  consider  an  arbitrarily  shaped 
flat  scatterer .  Take  N  rectangular  segments  to  approximate 


Fig.  4  A  flat  conducting  scatterer  approximated 


by  rectangular  subsections 


Add  (Ne-N)  segments  (dummy  segments)  to  get  a  full 
rectangle.  Now,  on  each  segment  (or  group  of  segments)  we 


choose  two  independent  current  expansion  functions;  one  in 


the  x  direction  and  the  other  in  the  y  direction,  as  shown 
in  Fig  .  5  . 


Fig.  5  Expansion  functions  used  on  the  n  subsection 


It  is  apparent  that  we  then  have  the  following  matrix 
equation . 


CZXX] 

[Zxy] 

Ix 

vx 

[Zyx] 

[Zyy] 

Iy 

vy 

(12) 


Here,  [ZXX],  [Z*y],  [ZyX],  and  [Zyy]  are  all  block  Toeplitz 

matrices  of  size  N  .  Renumbering  the  segments  in  terms  of 

e 

rows  and  columns,  (12)  can  be  rewritten  as 


Hj  »,  N?  », 

Ey  zxx  jX  +  V*  y  zxy  ^y 

m.a.n.n,  in,n_'*'  m,  nun,  n,,  n,n0 

n  -1  n  “1  1  2  1  2  1  2  n  "1  n  =1  1212  12 


m1  m2 


N2  n, 


m2  m, 


y  y  zy*  tx  +y  v  zyy  j 

m  in.n.n,  n  n_  *— •  m1m_n1n0  n.n, 

tlt«l  fl,»l  12  12  12  y,.,  12  12  12 


mim2 


m  si  ,2,. . .  ,N-| ;  m2  =  l  ,2  , . . .  ,M. 


(13) 


where  N1N2=Me,  N.,  being  the  number  of  segments  in  the  x 
direction  and  N2  being  the  number  of  segments  in  the  y 


direction.  Also 


,xx  _xx ,  . 

Z_  =  Z  ( m .  -n.  ,m0-n0  ) 

m  m 2n]n2  1  12  2 


!y  =Z  y  ( m  -n.  ,m0-n.  ) 

m1"2"ln2  '  1  2  2 


2S, 


Zyy  =zyy( m, -n. ,m0-n9 ) 

m1 m2n1 n2  1  1 ’  2  2 

Therefore  (13)  can  be  rewritten  as 


h2  M, 

£  2  Zxx(m1-n1 ,m2-n2)IX(n1 ,n2)+ 

n2  =  1  n  1  =  1 


h2  H, 


(14) 


y  Zxy(m1-n1  ,m2-n2>Iy(  n1  ,n2)  =VX(m1tm2) 


n2  =  1  n i  =  1 


i 


N2  ■  , 


T  Zyx(mrn1,n2-n2)IX(n1.n2)  + 


n2=1  n 1 =1 


n2  Mi 


Z2  Zyy(m1 -n1 ,m?-n2 ) I y( n1 ,n2 )  =Vy(m1tm2) 


r»2  =  1  n  i  =  1 


m1  si  , 2 , . . . , N1  ;  m2  =  1  ,2 , . . . ,H2 


(15) 


Now  (15)  is  obviously  a  convolution  equation  set,  written 
symbolically  as 


Z  *  I  +  Z  y  »  Iy  s  v 


Zyx  *  Ix  ♦  Zyy  »  Iy  =  Vy 


(16) 


For  the  second  example  consider  an  arbitrarily  shaped 
solid  imperfect  conductor  or  dielectric.  Now  take  N 
rectangular  cubic  segments  to  approximate  its  shape,  as 
shown  in  Fig  .  6 . 


Fig.  6  A  solid  dielectric  or  imperfectly  conducting 
scatterer  approximated  by  cubic  subsections 

Add  N@-N  dummy  segments  to  get  a  full  rectangular  cube.  To 
solve  the  problem  of  scattering  from  the  imperfect  conductor 
or  dielectric  using  the  MOM  formulation,  we  choose  for  each 
segment  (or  group  of  segments)  three  independent  current 
expansion  functions;  one  in  the  x  direction,  another  in  the 
y  direction  and  the  third  in  the  z  direction  as  shown  in 


segment 


Fig.  7  Expansion  functions  used  in  the  n  subsection 

It  is  apparent  that  we  then  have  the  following  matrix 
equat ion . 


>1 


[Zxx] 

[zxy] 

tzxz] 

[ZyX] 

[zyy] 

[Zyz] 

[Zzx] 

[Zzy] 

Czzz] 

[Zxx] 

.  [Zxy] 

,  etc 

(17) 


matrices  of  size  Ng.  Renumbering  the  segments  in  terms  of 
rows  and  columns,  we  can  write  (17)  as 


n3  n2  n1 

e  e  £ z™; 


XX  _  X 

„  I  + 

m  m^m^n^ngn^  n^n^n. 


nr 


m 


•a 


I 


O 


Z_>  L-*  L-t  n.)n2n3 

n3=l  n2=1  n 1 = 1 

n3  n2  n1 

EV  V*  zxz  iz  =  vx 

21...J  j  ^  1  ^2  ^2^3  n 1 n2n3  w ^  tn^  in* 

V'  n.,.1  n,=1 
etc  . 


(18) 


Here  N1N2N3=Ne,  N1  are  the  number  of  segments  in  the 
direction,  N2  are  the  number  of  segments  in  the  y  direction 
and  N3  are  the  number  of  segments  in  the  z  direction.  Also 


^ ^  tn^rn^  n  ^  n ^ ^ 3 


x  x 

=  Z  (m1 -n1 ,m2-n2 ,m3-n3 ) 


i1m2m3n1n2n3  =  Z*  U1  "n1  •  m2~n2  » m3”n3 ) 


m^  m2  m3  n  ^  n2  n3 


etc . 


—  Z  (w^ -n^  ,m2~n2 ,  m  3  ™  n  3  ) 


Therefore  (18)  can  be  rewritten  as. 


M3  h2  h, 


Y2  53  JZ  Zxx(m1-n1  ,m2-n2 ,m3-n3)  Ix(n1tn2,n3)  + 


n3=1  n2=l  n1=1 

N3  M2  N1 


7^  7~*  zyX(m1  -n1  .m2-n?  ,m  -n  )  Iy(nltn?,n  ) 


n3*1  n2  =  1  n  1  =  1 


15 


r>3  =  i  n2  =  1  n  1  =  1 

=  VX(m1 ,m2 ,m3)  (20) 

etc  . 

Now  (20)  is  obiously  a  set  of  convolution  equations. 
Symbolically,  we  can  write  this  set  as 

■Jxx  *  "^x  +  "£xy  *  ~ y  +  -|xz  v  "|z  _  ~ x 

Zyx  *  "ix  +  zyy  «"ly+*Zyz»"l2  =  Vy  (21) 

Zzx  *  1X  +  Tzy  *  ly  +  “z2Z  »  l2  =  Vz 

Here  denotes  three  dimensional  convolution. 

Other  examples  of  "two  and  three  dimensional"  problems 
of  this  type  (giving  sets  of  two  and  three  dimensional 
convolution  equations)  are  apertures  in  an  infinite  flat 
conductor,  antenna  arrays  with  both  polarizations,  non 
planar  antenna  arrays  with  more  than  one  polarization,  etc. 
Also,  some  problems  of  type  (1)  can  be  reformulated  as  type 
(2)  problems.  For  example,  a  planar  array  problem  using 
more  than  one  expansion  per  antenna  (say  two)  could  be 
written  in  the  form  of  (13)  and  hence  (16),  by  replacing  x 
and  y  in  the  equations  by  1  and  2  and  by  numbering  the  two 
expansions  1  and  2. 
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III.  SOLUTION  METHODS 

The  methods  that  we  use  to  solve  (6)  for  the  problems 
of  the  first  type,  and  (16),(21)  etc.  for  the  problems  of 
the  second  type  are  iterative.  We  will  discuss  here  in 

detail,  the  method  for  solving  (6).  To  demonstrate  the 
general  approach  for  the  second  type  of  problem,  we  also 

discuss  the  method  for  solving  (16). 

In  (6),  we  are  taking  the  convolution  product  of  the 
left  hand  side  and  equating  it  to  the  right  hand  side  for 

the  values  of  p^pg^.p^  in  the  re8ion  Sg .  However,  the  full 

convolution  process  given  by  the  left  hand  side  of  (6) 
produces  results  not  only  for  Sg  ,  but  also  for  the  regions 
outside  of  S_.  Specifically,  convolution  results  are 

C 

produced  for 

P1  =  -H1*2,-H1+3 . -1,0, 1,2 . 2N1-1 

P2  S  -n2+2 »-n2  +  3 » • • • •”1  »° *  1  >  2 » • • • , 2N2-1  (22) 

P3  =  -N  3+2 , -N  3+3 . -1  ,0,1  ,2,. . . ,2N3-1 

etc  . 

However,  values  of  V  not  in  the  region  S  are  unknown  and 

values  of  J  not  in  the  region  S  are  known  (equal  to  zero). 

If  we  call  the  values  of  "v  in  region  S  to  be  "v*  (for 
impressed)  and  outside  S  to  be  V  ,  then  (6)  can  be  rewritten 
as 

«■*  i  *4Q 

Z»J  =V>V  =  V  (23) 

Here,  no  restrictions  are  placed  on  the  region  of  validity. 


If  we  take  the  Discrete  Fourier  Transform  (DFT)  of 


(23),  on  the  basis  of  3N1-2  elements  for  p1 ,  3N2-2  elements 
for  p2  etc.,  then  we  get  an  algebraic  equation  [2]. 

Z  J  =  V  (24) 

Here  "  ***  "  denotes  transformed  quantities.  Equation  (24)  is 
true  for  each  transformed  quantity;  in  other  words, 

^('<^.1<2,***,*<m^^1,^2 . =  ^ 

k  =1,2 . 3N  1  -2  ;  k?  =  1  ,2 . 3M2-2  ;  .  .  .  ;  etc.  (25) 

Therefore,  if  we  know  V  for  all  values  of  p  p2>..pM,  then  we 
can  determine  V  and  find  J  by 


V(k. ,k2 , . . . ,kM) 


J  ( ,  k 2  ,  *  «  •  » k^ ) 


(26) 


Z(  k1  ,  k  2  ,  •  t  »  *^Jy|) 

The  inverse  DFT  ( IDFT )  then  gives  7.  However,  we  know  only 


V1  and  not  V°, 


Thus  the  following  procedure  is  used; 


STEP  1-Assume  V 


Normally  we  take  all  (initial)  values 


of  *$°  to  be  zero.  (As  shown  later  in  the  Appendix,  the 
"distance"  of  the  initial  guess  from  the  correct  value 
does  not  effect  the  convergence;  only  the  number  of 
iterations  needed.)  Call  this  first  guess  of  V°  by 
STEP  2-Take  the  DFT  of  Z  on  the  basis  of  3^-2  for  p1  , 
3N 2~2  for  p2,  etc.  to  get  Z. 


STEP  3-Compute  V 


(1  ) 


**  ♦*?,> 


(27) 

STEP  4-Take  the  DFT  of  V(1)  on  the  same  basis  as  in  step 


2  to  get  V(1)(k,  ,kj . kM> 

STEP  5-Compute  7^  1 ^ ( k1 , k2 » . • . . kM)  using  (26) 


* 


1iV 


m 


s 


vi 


STEP  6-Take  the  IDFT  (on  the  same  basis  as  DFT  in  step  2) 
of  J(1)  to  get  J ( 1 ) . 

STEP  7-Since  is  not  the  correct  answer,  it  will  have 

nonzero  values  outside  S.  Change  the  values  of  ^ 
outside  S  to  zero.  (This  is  the  same  as  truncating  or 

*■*  -•  n 

projecting  onto  S. )  Call  this 

STEP  8-Take  the  DFT  of  )  to  get  jj>1)(k1,k2 . kM) 

STEP  9-Compute  (  =  Z  JP,  as  given  in  equation  (25)). 


STEP  10-Take  IDFT  of  V1 
VP=Z»JP  .) 


to  get  V 


(2)  • 


(Note 


STEP  11-Since  ^(-|)  is  n°t  yet  the  correct  answer,  values 

d  ^  i 

of  on  S  are  not  be  equal  to  V  .  Here,  we  can  check 

the  accuracy  by  comparing  "v1  with  on  S.  One  method 

is  to  check  the  maximum  ("v  i-V  ^  )/"v 1  for  all  elements,  as 

well  as  the  average.  If  the  maximum  and  average  are 

below  a  certain  value  (say  .1%  and  .01%  respectively), 

then  stop.  We  can  also  use  the  criterion  of  the 

convergence  of  ^(n)»  i.«..  the  relative  magnitude  of 

J(n)-J(n  i)  in  comparison  to  J(n)*  or  combine  the  two 

criterion  s . 

STEP  12-If  we  decide  that  more  iterations  are  needed, 

n  i 

then  change  the  values  of  on  S  to  V  .  Call  this 

7(2)‘ 

STEP  13-Replace  in  step  4  by  V(2)* 

STEP  14-Continue  from  step  4  onwards  until  the  criterions 
in  step  11  are  satisfied. 


W5? 


The  solution  techniques  for  sets  of  convolution 
equations  produced  by  the  problems  of  the  second  type  are 
not  as  straight  forward.  To  illustrate  the  general 

technique,  consider  (16).  In  (16),  we  are  taking  the 
convolution  products  of  the  left  hand  sides  and  equating 
them  to  the  right  hand  sides  for  values  of  in  the 

region  S@.  However,  the  full  convolution  process,  as  given 
by  the  left  hand  sides  of  (16)  produces  results  not  only  for 
S_ ,  but  also  for  regions  outside  S. .  Specifically, 

convolution  results  are  produced  for 
m1  =-N 1 +2 , -N 1 +3 . . • . . -1  .0,1,2 . 2N1-1 

n»2  =-N2+2,-N2+3 . -1,0, 1,2 . 2N2-1  (28) 

^  X 

as  in  the  first  type  of  problems.  However,  values  of  V  and 
Vy  not  in  the  region  S  are  unknown,  and  the  values  of  "lx  and 
Iy  not  in  the  region  S  are  known  to  be  equal  to  zero.  If  we 
denote  the  values  of  Vx  and  Vy  in  the  region  S  by  Vx*  and 
Vyi  (for  impressed),  and  those  outside  S  by  Vxo  and  Vyo, 
then  (16)  can  be  rewritten  as, 

ZXX»IX  +  Zxy*Iy  =Vxi+Vxo  =VX 

Zyx»lx  ♦  Zyy»Iy  =Vyi+Vyo  =Vy  (29) 

Here,  no  restrictions  are  placed  on  the  region  of  validity. 

If  we  take  the  DFT  of  (29)  on  the  basis  of  3N 1 —2 
elements  for  m1  and  3N2-2  elements  for  m2 , 
algebraic  equations  [2] 


we  then  get  the 


*\  v  V' 


'  V*  L"V 


2xx<k,.k2)rx(k1.kj).ry(k1,k2)iy(k1.k2).v,t<k,,i,2) 

?,x(k,.k2)ix(k1,k2).iyy(k1.k2>-iy(k1,k2)4''(k1.k2) 


^  1  =  1  *  2  ,  . . . , 3N  ^ -2 ;  k2  =  1  , 2 , . . . , 3N  2— 2 


(30) 


Equation  (30)  is  true  for  each  transformed  quantity. 

Therefore,  if  we  know  Vx  and  Vy  for  all  values  of  m  and  ra 

1  2  • 

we  can  find  V  and  Vy.  Since  (30)  can  be  written  as. 


zxx(k1,k2)  zxy(kvk2)  px(k1,k2)  |vx(k1,k2) 


?,X(k,.k2)  VH  k,,k2)  iy(k,,k2)  vy(k,.k2) 


(31) 


we  can  find  lx  and  ly  easily  as 


w  x  , 

I  (k1 ,k2) 


Zyy(ki  ,k2)vx(k1  ,k2)-zxy(k1  ,k2)vy(k1  ,k2) 

Zyy(k1  ,k2)lXX(k1  .k^-Z^U,  ,k2)ZyX(kl  ,k2) 


(32) 


ly<k1fk2)  = 


vx(k1fk2)  -  zxx(k1,k2)Tx(k1,k2) 


zAy(k1fk2) 


The  IDFT  then  gives  IX  and  "fy.  However,  we  know  only  Vxi, 
Vy  an^  not  V  ,  Vy  .  Thus  the  following  iterative 
procedure  can  be  used. 

STEP  1-Assume  Vxo,  \7yo  .  Normally,  we  take  all  elements 


of  V  ,  VJ  to  be  zero.  Denote  this  first  guess  of  V  , 
Vyo  by  VXO  V° 

STEP  2-Take  the  DFT  of  1XX,  Zxy,  ZyX,  and  lyy  on  the 
basis  of  3N 1 -2  for  m1  and  3N2-2  for  m,  to  get  ZXX,  7xy, 
Zy\  and  *Zyy. 


-•x  *-•  x  i  -•  xo 

STEP  3-Compute  Vd)  =  v  +  v(  i : 


(33) 


;yi  ♦  vy°. 


STEP  4-Take  the  DFT  of  V 


(1  )  ’  v(  1  ) 


on  the  same  basis  as  in 


step  2  to  get  Vx1j(k1,k2)  and  Vy ^ ^ ( k1 ,k2) . 

STEP  5-Compute  lx 1  )  ( k,  ,  k2 )  ,  I^1)(k1,k2)  using  (32). 

STEP  6-Take  the  IDFT  (on  the  same  basis  as  DFT  in  step  2) 
of  IXn  and  lf1}. 

STEP  7-Since  "T^  j  and  j  ,  are  not  the  correct  answer, 
they  have  nonzero  values  outside  S.  Therefore,  change 


the  values  of  I, 


and  I; 


outside  S  to  zero.  (This  is 


U II C  V  dl  UtJ  3  X  £  ^  J  dliu  ^  ^  u  ^  i  w  •  v  i  ii*a 

the  same  as  projecting  7X 1  ^  and  ^  onto  S.)  Call  these 

- x p  .  -yp 

X(1)  and  X(  1  )  • 

STEP  8-Take  DFT 1  s  of  IXf).  to  get  I*f)(k1tk2)  and 

TyP  ( k  v  ' 

(1)1*2' * 

STEP  9-Compute  Vx£j  and  as  given  in  (30). 

STEP  10-Take  IDFTs  of  V*2)  and  7]^  to  get  'vX|) 


and  7^..  (Note  that 


V*P  /7ltx.T1(P  +lxy»1yp  and 
V(2)"Z  I(2)  +  Z  X(2)  anG 

«yp  -lyx»ixp  +'zyy*Typ  ) 

V(2)-Z  1 ( 2 )+Z  X(2)  ' 

STEP  11-Since  I^  ^  ,  J(i)  are  not  Yet  the  correct  answers, 
values  of  Vx|)  and  "vy|)  on  S  are  not  equal  to  dxi  and 
Vyi.  Here  we  can  check  the  accuracy  by  comparing  “Vx*, 


wi  th 


on  S.  The  same  kind(s)  of 


criterion(s)  as  in  step  11  of  the  solution  procedure  for 
the  problems  of  the  first  type  can  be  used  to  determine 
whether  or  not  the  iteration  has  converged, 

STEP  12-If  we  decide  that  more  iterations  are  needed, 
then  change  the  values  of  VXp^  and  Vyp^  on  S  to  be  Vxi 


and  V 


Denote  these  V*,),  Vy?j. 


STEP  IB-Replaced*^  andd^^  in  step  4  by~V*?)  and 


STEP  14-Continue  from  step 


onwards  until  the 


criterionC s)  of  convergence  in  step  11  are  satisfied. 


IV.  SAMPLE  COMPUTATIONS  AND  COMMENTS 


Computer  Programs  using  the  techniques  devised  in  the 
preceding  sections  have  been  written.  They  are  listed  in 
the  Appendix  of  this  report.  In  this  section  we  give  the 
results  of  computations  that  were  made  using  these  programs. 
The  routines  and  formulations  for  computing  the  impedence 
matrix  (or  mutual  coupling  matr  i  x )  (z)  are  from  [33.  C4],  and 
[53. 


Fig.  8.  Scattering  from  a  thin  straight  wire. 

Table  1.  Comparison  of  DCM  and  LU  decomposition  methods 


L 

1 

N 

s 

Ne 

Nc 

LU  decomp 

(mult,  ops  . 

I 

) 

DCM 

(mult . 

0.75 

0.0625 

1  2 

5 

1  5 

42 

10 

12  80 

2.0 

0.05 

40 

1  9 

64 

2286 

9 

6912 

4.0 

0 . 1 

40 

1  9 

64 

2286 

6 

4608 

8.0 

0.1 

80 

39 

128 

19773 

1 1 

19712 

16.0 

0.1 

160 

79 

256 

1 64346 

5 

20480 

23.85 

0.075 

318 

158 

512 

1314771 

12 

135168 

32 .0 

0.1 

320 

159 

512 

1339893 

6 

67584 

Here  a  =  wire  radius  in  wavel engthsC 0 . 01 3477089  for  problems 


in  Table  1 ) 


L  =  length  of  wire  in  wavelengths 

1  =  length  of  each  segment  in  wavelengths 

Ns=  number  of  segments 

Ng=  number  of  expansion  functions  needed 

N  =  basis  on  which  FFT  and  IFFT  are  taken  for  DCM 
c 

I  =  number  of  iterations  needed  to  get  maximum  error  in 
current  to  be  less  than  1%  or  maximum  error  in 
field  to  be  lc33  '..hail  .  1? 

From  Table  1  we  see  that,  based  on  the  number  of  complex 
multiplications  required,  DCM  is  faster  than  LU  decomposition 
for  problems  with  more  tuan  AO  expansions.  Total  computing 
time  needed  to  solve  the  32  wavelengths  problem  Clast  entry 
in  Table  1)  is  14.16  seconds  on  an  IBM  4341.  This  includes 
computing  time  needed  to  set  up  the  impedance  matrix. 

The  problem  of  two  thin  wire  scatters  with  a  gap  in 
between,  as  in  Figure  1,  was  also  solved.  The  number  of 
iterations  needed  for  .123%  maximum  field  error  and  .00404% 
average  field  error  was  found  to  be  14.  The  problem  is  the 
same  as  the  23.85  wavelengths  wire  problem  of  Table  1, 
except  that  114  segments  in  the  middle  are  missing.  Since 
the  original  problem  needed  12  iterations  for  the  same  level 
of  accuracy,  the  insertion  of  the  gap  does  not  seem  to  cause 
much  increase  in  computing  time. 


Fig.  2  shows  the  problem  of  radiation  from  a  linear 
antenna  array.  If  one  expansion  function  per  antenna  is 
used,  then  the  Method  of  Moments  formulation  gives  a  matrix 
equation  equivalent  to  a  one  dimensional  convolution 
equation.  The  formulation  used  is  as  given  in  [4],  Since 
the  mutual  coupling  matrix  is  Toeplitz  in  this  case,  it  can 

p 

be  solved  using  the  faster  (N  order)  algorithm  for  Toeplitz 

matrices  as  given  in  [43.  Therefore,  Table  2  compares 

between  computing  time  needed  for  DCM  with  the  computing 

2 

time  needed  for  the  ( N  order)  Toeplitz  algorithm  given  in 
[43.  The  computing  time  measurements  were  made  on  an  IBM 
4341.  The  matrix  set  up  time  is  not  included,  which  would 
be  the  same  for  both  cases.  All  the  problems  in  Table  2  are 
linear  arrays  with  halfwave  antennas  one-quarter  wavelength 
in  front  of  an  infinite  ground  plane.  The  seperation 
between  antennas  is  one-half  wavelength  also.  Uniform 
excitation  is  used  for  all  cases  given  in  Table  2.  Other 
excitations  were  tried  and  number  of  iterations  needed  (and 
hence  computing  times)  were  found  to  be  practically 
independent  of  the  type  of  excitation. 

We  see  that,  for  very  large  arrays,  DCM  is  faster.  For 
1000  antenna  elements  it  is  nearly  4  times  faster. 
Breakeven  seems  to  occur  at  about  300  antenna  elements. 
Comparisons  between  the  solutions  given  by  DCM  and  Toeplitz 
algorithms  were  made  for  all  problems  except  the  1000 
antennas  problem,  and  the  agreement  was  to  within  .11 
maximum  difference  in  current  in  all  cases. 
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Table  2.  Comparison  of  DCM  and  Toeplitz  methods 


"e 

I 

Setup 

DCM 

Toeplitz 

Field  Error 

Time 

Time 

Time 

(  sec) 

(  sec ) 

(  sec ) 

(%) 

1 1 

4 

.57 

.37 

44 

3 

1.55 

1.23 

88 

3 

1  .51 

2.43 

1  .57 

.0222 

. 000971 

250 

3 

3.97 

9.8 

9.79 

.0224 

.000349 

352 

3 

5.55 

19.48 

18.53 

.022 

.000935 

1000 

3 

11.57 

34.22 

135.44 

.0218 

.000423 

1000 

4 

- 

45.63 

- 

.00466 

.000365 

Here  N 

e 

is 

the  number  of 

antennas  in 

the  array 

I  is  the  number  of  iterations  needed  for  the  given 
field  error  and  the  last  current  change.  For  both  field 
error  and  last  current  change,  the  upper  entry  is 
maximum  error  and  the  lower  entry  is  average  error. 


N 

-  V  ^ 


Vs-' 
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Computing  time  savings  are  even  more  dramatic  for 
antenna  arrays  with  some  antennas  missing,  i.e.,  gaps.  The 
matrix  produced  by  MOM  in  this  case  is  no  longer  Toeplitz, 
and  LU  decomposition  (needing  1/3  N  ^  multiplicative 
operations)  is  usually  used  instead  of  Toeplitz  methods 
(needing  2N  multiplicative  operations).  But,  as  explained 
in  section  II,  we  can  add  dummy  segments  and  3till  use  DCM 
to  solve  the  problem.  The  number  of  iterations  needed 
increased  by  only  1  over  that  needed  for  the  same  problem 
without  gaps  in  each  of  the  cases  that  were  tried.  The 
results  are  given  in  Table  3. 

Table  3.  Results  for  array  problems  with  gaps 


N  N 

e  g 

Gap(  s) 

I 

Field  error(t) 

44  1 

17-28 

4 

.0166 

.  00 1  75  (  av  g) 

2 

12-16 

4 

.0384 

28-32 

. 00759 ( avg) 

3 

12-15 

4 

.0330 (max ) 

21-26 

.00542 

33-37 

Here  Ng  is 

the  number 

of  gaps 

Gap(  s) 

gives  the 

start  and 

end  segment  numbers 

gaps 

I  is  the  number  of  iterations  needed  for  the  given 


field  error 
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With  the  LU  decomposition  method,  the  computing  time  would 
be  10  times  larger  even  for  the  *14  antenna  problem. 

For  linear  antenna  arrays  which  are  lined  up  at  an 
angle  as  shown  in  Fig.  9,  then  one  expansion  per  antenna 
would  no  longer  be  enough  since  the  current  will  not  be 
symmetric  anymore. 


Fij’.  9  Linear  antenna  array  lined  up  at  an  angle 

The  MOM  formulation  will  no  longer  give  an  impedance  matrix 
that  is  Toeplitz  but  will  give  an  impedance  matrix  that  is 
block  Toeplitz.  This  is  equivalent  to  the  two  dimensional 
discrete  convolution.  We  can  solve  by  using  the  two 
dimensional  DCM  technique  but  since  one  of  the  dimensions 
will  have  only  three  points,  it  is  not  worthwhile.  However, 
we  can  look  at  the  matrix  equation  as  three  one  dimensional 
convolution  equations  similar  to  what  is  done  in  (17)  to 
(21).  However,  convolutions  would  be  one  dimensional  here, 
instead  of  three  dimensional  as  in  (21).  If  we  call  the 


first,  second  and  the  third  expansions  of  each  antenna,  a,b, 
and  c  respectively,  then  the  equivalent  set  of  equations  is 


zaa 

ft 

la 

+ 

ab 

> 

Ib 

-4 

zac 

• 

G 

IG 

-* 

=  0 

zba 

ft 

Ia 

4 

zbb 

ft 

Ib 

4 

Zbc 

• 

Ic 

=  V 

(3M  ) 

zca 

ft 

Ta 

4 

zcb 

• 

Tb 

4 

zoc 

Tc 

-• 

=  0 

The  solution  of  (34)  using  DCM  would  require  nine  times 
more  multiplicative  operations  than  the  solution  of  (9). 
However,  since  the  block  Toeplitz  solution  technique  will 
also  need  nine  times  more  than  the  Toeplitz  solution 
technique,  the  timing  comparisons  of  Table  2  will  remain 
unchanged.  For  problems  with  gaps,  if  three  expansions  per 


antenna  are  used,  DCM  will 

still 

need 

only  nine 

times 

more 

multiplicative 

operation s 

per 

iteration  . 

But 

LU 

decom  po  si tion , 

being  1/3  N3, 

will 

need 

twenty 

seven 

times 

more  multiplicative  operations. 

Fig.  3  shows  the  problem  of  radiation  from  a  planar 
antenna  array.  If  only  one  expansion  function  per  antenna 
is  used,  then  the  MOM  formulation  gives  a  matrix  equation 
which  is  block  Toeplitz.  This  is  equivalent  to  a  two 
dimensional  convolution  equation.  The  formulation  is  as 
given  in  C 5 ] .  Since  the  matrix  equation  is  block  Toeplitz, 
it  can  be  solved  using  faster  (N  order  where  a  e2.5) 
algorithm  for  block  Toeplitz  matrices  as  given  in  [53. 
Table  4  lists  the  computing  time  requirements  for  the  DCM 
solution  of  some  planar  array  problems.  The  computing  times 


given  include  the  setup  times.  All  the  problems  are  for 


planar 

array  problems 

with  halfwave 

antennas  one-quarter 

wavelength  in  front 

of 

an  infinite 

ground  plane.  The 

seperation  between  antennas  is  one-half 

wavelength 

in  either 

direction. 

Table  4 

.  Results  for 

some 

planar  array 

problem  s 

Ne 

Excitation 

I 

Com  put i ng 

Field  Error 

Current 

Time( sec ) 

(X) 

Change 

(%) 

1  6 

Uniform 

4 

2 

.05 

.25 

(4x4  ) 

.03 

.14 

36 

Uniform 

4 

3 

.056 

.25 

(6x6) 

.025 

.1 

Exponential 

4 

3 

.055 

.27 

Taper 

.021 

.087 

Beam  Steer 

4 

3 

.034 

.62 

(45  *-45  • ) 

.008 

.056 

Beam  Steer 

4 

3 

.02 

.18 

(30--20 •) 

.007 

.05 

Progressive 

4 

3 

.074 

.3 

Phase  shift 

.032 

.  1 

( 30  * -20  ’ ) 

121 

Uniform 

4 

13 

.07 

.29 

.03 

.08 

Exponential 

4 

5 

.068 

.24 

N  Excitation 


18*4  9  *  Uniform 


Beam  Steer 


I  Computing  Field  Error  Current 


Time (sec ) 


Change (% ) 


<30*— 30”) 


.009 


.0014 


Here  is  the  number  of  antennas  in  the  array 

I  is  the  number  of  iterations  needed  to  get  the  given 
accuracy.  For  both  field  error  and  (last)  current 
change,  the  upper  entry  is  the  maximum  and  the  lower 
entry  is  the  average. 


He  see  that  even  for  very  large  arrays  DCM  is  quite 
fast.  The  problem  with  1849  antennas  takes  only  4  minutes 
computing  time.  Just  as  for  the  linear  arrays,  problems 
with  gaps  could  also  be  treated  and  would  take  comparable 
amounts  of  time. 


For  more  accurate  planar  array  solutions,  three 
expansions  per  array  should  be  used  since  the  current  on 
each  antenna  is  not  symmetric.  This  would  increase  the 
required  computing  time  by  a  factor  of  nine  as  before. 
Notice  also,  the  independence  of  the  number  of  iterations 
required  to  the  size  of  the  array. 
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To  check  the  accuracy  of  the  DCM  technique,  the 
solution  of  the  36  antenna  planar  array  was  made  using 
matrix  inversion  routine  LINEQ  from  [3],  and  the  agreement 
between  the  currents  were  found  to  be  better  than  .IX 
maximum  difference. 

The  problem  of  a  planar  array  with  antennas  arranged  in 
diamond  patterns  instead  of  rectangular,  can  also  be 
formulated  as  a  two  dimensional  convolution  equation  by 
adding  dummy  elements  (as  shown  in  Figure  11),  to  make  a 
parallelogram.  The  diamond  pattern  arrangement  is  shown  in 
Fig.  10. 


Fig.  10  The  Diamond  pattern  arrangement 


Fig.  11  A  planar  array  with  antennas  arranged  in  diagonal  patterns 
The  MOM  formulation  using  one  expansion  per  antenna 


will  then  give  a  block  Toeplitz  matrix  which  can  be  solved 
using  two  dimensional  DCM.  However,  it  cannot  be  solved 
using  the  block  Toeplitz  method  since  the  field  on  the  dummy 


elements  are  unknown. 

Therefore,  only 

LU 

decompo  si  tion 

or 

two 

dimensional  DCM 

can  be 

used  . 

For 

a  large  array. 

DCM 

will 

be  considerably 

faster  . 

Since 

the 

current  on 

each 

antenna  is  not  symmetric,  using  three  expansion  functions 


per  antenna  will  give  a  much  more  accurate  result  and  will 
need  nine  times  more  computing  time  for  the  DCM.  With  LU 


decomposition  method  computing  time  will  go  up  twenty  seven 


times  the  already  large  value. 

V.  DISCUSSIONS 

Discrete  Convolution  Method  for  solving  the  matrix 
equation  set  up  by  the  MOM,  is  found  to  be  accurate  and  much 
faster  than  either  the  Gaussian  Elimination  (LU 
decomposition)  or  the  Toeplitz  and  block  Toeplitz  methods 
given  in  [4]  and  [5].  Also  DCM  can  solve  a  wider  range  of 
problems  than  the  Toeplitz  and  the  block  Toeplitz  methods. 
The  number  of  iterations  needed  by  DCM  for  a  given  accuracy 
is  also  found  to  be  practically  independent  of  size. 
However,  it  is  dependent  on  other  factors.  For  instance, 
the  number  of  iterations  needed  to  solve  an  array  problem  is 
found  to  be  dependent  on  antenna  length,  antenna  seperation 
and  the  ground  plane  distance. 

With  careful  formulation,  the  problem  of  radiation  from 
a  planar  array  backed  by  a  finite  ground  plane  can  be  solved 
with  the  DCM.  Therefore  DCM  may  prove  to  be  useful  in 
designing  array  antennas. 

The  other  numerical  technique  using  FFT  and  IFFT  to 
solve  iteratively,  electrically  large  problems  is  the 
Spectral  Theory  of  Diffraction  (STD).  This  technique  is 
well  known  and  a  number  of  papers  and  reports  [6],  E73.C8] 
etc,,  have  been  published  about  STD.  The  difference  between 


STD  and  DCM  is  that  STD  solves  the  problems  in  the  spectral 
domain  and  DCM  solves  the  problems  in  the  spatial  domain. 
Therefore,  for  certain  types  of  problems  STD  may  feel  more 
natural  and  for  other  types  of  problems  (for  example  planar 
arrays),  DCM  may  feel  more  natural.  Also  errors  in  each 
technique  have  different  causes. 

Since  both  STD  and  DCM  are  iterative  techniques  using 
FFT  and  IFFT,  they  both  have  numerical  inaccuracies 
associated  with 

(i)  the  use  of  FFT  and  IFFT 
(ii)  taking  only  a  finite  number  of  iterations  based 
on  some  criterion 

However,  as  shown  in  the  Appendix,  for  DCM  these  errors  are 
insignificant.  But  since  DCM  is  the  iterative  solution  of 
the  MOM  formulation  of  the  original  problem,  DCM  will  have 
in  addition  the  inaccuracies  associated  with  the  MOM 
formulation  (but  not  the  matrix  inversion).  On  the  other 
hand,  STD  has  the  following  numerical  problems  [6], 

(  i )  wi ndo wi ng  and 

(ii)  the  need  to  take  sufficient  number  of  points  to  make 
sure  that  the  aliasing  effect  is  small 

Thus  the  numerical  errors  of  DCM  and  STD  are  of 
different  natures.  However,  since  the  MOM  formulation  has 
been  in  wide  use  for  a  considerable  period  of  time,  the 
numerical  errors  associated  with  the  MOM  are  familiar 
through  experience. 


APPENDIX 


I.  INDEPENDENCE  OF  CONVERGENCE  ON  STARTING  POINT 


We  prove  here, 

the 

independence 

of  convergence 

on 

starting  point  for 

the 

one 

dimensional  case. 

Proofs 

for 

higher  dimensions 

will 

be 

similar . 

Consider 

the 

one 

dimensional  problem 

given 

i  n 

(9).  The 

discrete 

convolution 

equation  to  be  solved  is, 

Z  *  J  =  V  ( A1  ) 

But  V  is  unknown,  only  B[V]  is  known.  Here  the 
function  8[  ]  means  truncate  values  outside  region  S  and 
replace  with  zeros.  We  also  know  that  J  is  confined  to  S, 
i  .e  . 

J  =  8  [ J ]  ( A2 ) 

The  general  solution  technique  is, 

Z  »  Jn+1  =  SCV]+  6  [?•  8  [Jn]]  ( A3 ) 

/•*  - 
where  6[  ]  is  the  complement  of  “[  ], 

-»  -*  th 

Jn  is  the  approximation  of  J  at  the  n  iteration 

The  correct  solution  is, 

Z  *  J  =  ©  C  V  3  4-  0  [V]  ( A4 ) 

Therefore  from  (A3)  and  (A4),  we  get 

Z  »  ( J-Jn+1 )  =  6  CV-Z  »  6  C  Jn3  3  ( A5 ) 

Using  (A1),  the  equation  above  becomes, 

l*  (J-J.  .)  =  8tZ»T  -  Z« 6  [J]]  ( A6 ) 


*£&&& 


l~.  J-J  W7WZ  f  ,-T' 


Kt-i 


KV 


rvr»."  *.  %  ‘.  V".  /.  ,-s  ,-» 
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But  by  ( A2) , 


z  »  ( j-Jnvl )  =  e [z*  6 [ j-j ni i 


(  A7) 


Therefore,  if  we  denote  the  error  in  the  approximate 
th 


solution  at  n 


step ,  (J-J  )  as  £ 
r  n  n 


z  »  £ 


n+1 


=  e cz*  e [£n]] 


(A8) 


We  now  prove  the  independence  of  convergence  to  the 
starting  point,  in  the  sense  that  if  the  convergence  is 
achieved  for  a  certain  starting  point,  then  the  convergence 
is  achieved  for  other  starting  points  nearer  or  further  than 
that.  Suppose  that  we  achieve  convergence  if  we  start  with 


the  error  £q.  If  instead  we  start  at  a  different  starting 


point  with  the  error,  SQ=  *  tQ ,  then 


Z  *  S. 


e  [z*  e  [<*£o3] 


z  «  S1  =  e  [Z*  6  [ £0n 


z  *  sl  =  *  z  *  e1 


( A9) 

( A10  ) 
(All) 


Therefore , 

-  7 

Z  •  (S1-*E1  ) 


( A12  ) 


From  the  fact  that  the  original  convolution  equation  is 
the  Method  of  Moments  formulation  of  the  physical  problem 
which  can  have  no  currents  for  zero  excitation,  (A12)  can  be 
interpreted  as  indicating, 


1 


1 


( A13) 


and 


*  £ 


(  A 1  4  ) 


►  T 


:: 


■f--*  VA.VV' 


Since  the  convergence  is  achieved  when  we  start  with  the 


error  £q,  in  the  limit  as  n  approaches  00  ,  fin  approaches 
zero.  Therefore, 

lim  $n  =  0  ( A 1 5 ) 

r\~*  to 


II.  CONDITION  FOR  CONVERGENCE 

The  condition  for  convergence  for  the  one  dimensional 
case  is  given  here.  Multi-dimensional  cases  will  have 
similar  conditions  for  convergence.  The  general  solution 
technique  as  given  by  (A3)  for  the  one  dimensional  case  is, 

Z  •  Jn  +  1  =  &  [V  3  +&[Z»  ©  C  J  n  ]  1  (  A 1 6  ) 

However,  it  is  apparent  that  since  the  convolution  can  be 
written  as  a  matrix  multiplication  operation, 

t  yi  7n>i  =  9CV]  ♦8cc^]ecjn]]  ( A 1 7  ) 

Here  [  ^  ]  is  the  circulant  matrix  produced  from  ~i  and  not 
the  same  as  [Z]. 

The  trunction  operator  0  (  ]  can  also  be  represented  as  [T], 
a  diagonal  matrix  with  I's  at  places  on  the  diagonal 
corresponding  to  region  S  and  zeros  elsewhere.  Similarly, 

A 

it  is  easy  to  see  that  the  operator  0[  ]  can  be  represented 
* 

as  [T],  a  diagonal  matrix  with  I's  at  places  on  the  diagonal 
corresponding  to  region  S  (i.e.  the  region  outside  S)  and 
zeros  elsewhere.  Therefore  C  A 1 7 )  can  be  rewritten  as, 

C  J']  Jn+1  =  [T ]  V  [T][  J][T]  jfn  (  At 8 ) 

7n+1  =  [  J]_1  CT 3  V  ♦  [  £]_1  CT 3 C  JOCT]  Jn  (  A 1 9 ) 


goastfsife 
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Since,  C  is  the  circulant  matrix  with  Z  as  its  rows,  it 


has  an  inverse  if  Z  has  a  DFT  [9].  Let 
C Q 3  =  C  ^r1  ET3 


[  R  3  =  [  flrV  CT3C  4/3  [  T] 


Therefore , 


Jn+1  =  [Q3  V  +  [R]  Jn 


(  A20) 
(  A21  ) 


(  A22  ) 


It  is  easy  to  see  from  (A2)  and  (A19)  that  the  exact 


solution  is  given  by 


J  =  C Q ]  V  +  [R]  J  (A23) 

Equation  (A22)  can  be  rewritten  as, 

^n+1  =  CQ3  T  +  [R][Q]  7  +  [R]2  Tn  (A24) 

Repeated  application  of  (A22)  gives 

Jn+1  =  (C I3+[R]+[R]+.  .  .) [Q]  V  +  C  R  3  ^  Jn  (A25) 
If  the  maximum  eigenvalue  of  [R],  ^  max  is  such  that 


(  A24 ) 


Then  [9 ] , 

CI]  +  [R]  +  [R]2>. . .  =  C I  —  R 1 ~ 1 

and 

lRin+1  )  =  0 

as  n  approaches  oc . 


( A26) 


( A27 ) 


( A28 ) 


There  fore 


C I  —  R  3 —  1  CQ]  V* 


[I-R3  Jn+1  =  t Q 3  V 


=  [Q3  V  ♦  C R 3 J , 


as  n  approaches  00. 

Therefore  from  (A23)  and  (A31), 


(A29) 
(A30) 
( A3 1  ) 


as  n  approaches  oo}  if 

1^1  <1 

max 


(A33> 


However,  the  usefullness  of  the  above  condition  is  limited 

A 

as  the  computation  of  ^  requires  fp  order  complex 

multiplications  and  so  will  take  longer  than  the  solution 

itself,  although  not  as  long  as  using  Gaussian  Elimination 

o 

which  requires  complex  multiplications  of  the  order  N  . 


III.  ESTIMATION  OF  NUMERICAL  ERRORS 


The  Discrete  Convolution  Method  (DCM)  is  an  iterative 
technique  for  solving  the  matrix  equation, 


[  Z]  J 


(A34) 


formulated  by  the  Method  of  Moments.  This  is  done  by 
looking  at  the  above  equation  as  a  convolution  equation. 


Z  *  J 


(A35) 


Now,  it  is  apparent  that  if  we  are  given  the  answer  for 
(A31*).  (say),  then  we  can  use  J_  to  take  the  matrix 

3  3 

multiplication  with  [Z]  and  get  T  ,  which  we  can  then  check 

3 

against  V  to  see  if  the  answer  given,  Jg  is  correct  or  not. 
Also,  instead  of  taking  FFT  and  IFFT,  if  we  actually 


convolve,  by  using  the  relationship 

N 

V  =  22  Z ,  ,  J 

m  £77  (m-n)  n 


(  A36) 


to  compute  \T  from  J  ,  then  the  computations  involved  in 

3  3 

using  (A35)  is  identical  to  that  involved  in  using  (A31*). 


Therefore,  if  the  matrix  equation  solution  of  (A31*) 
using  Gaussian  Elimination  is  unique  (in  the  sense  that  up 
to  the  desired  precision  point  there  are  no  two  answers  to 
A31*,  although  there  may  be  many  beyond  that  precision 
point),  then  trying  out  various  J  ' s  (chosen  randomly,  found 
by  iteration  or  any  method)  in  (A31*),  would  also  give  us  the 
same  unique  answer.  Hence,  trying  out  various  J  's  in 

3 

(A35),  provided  we  actually  convolve  (i.e.  use  A36),  will 
give  the  same  answer.  So,  any  difference  in  DCM  and  matrix 
inversion  solutions  will  come  only  through 

(a)  the  use  of  FFT  and  IFFT 

(b)  the  fact  that  the  iteration  is  carried  out  only 
up  to  a  certain  accuracy  based  on  some  criterion. 


E(  error)  =  J~n/3 


Numerical  errors  due  to  (a)  can  be  analyzed  in  the 
following  way.  By  [2],  the  use  of  FFT  introduces  the  output 
error,  the  expected  value  of  which  is 

(A37) 

where  b  is  the  number  of  machine  precision  bits 
N  is  the  basis  on  which  FFT  is  made  and 
E(error)  is  the  expected  value  of  normalized  error 
Since  IFFT  needs  identical  computations  as  FFT,  IFFT  causes 
error,  the  expected  value  of  which  is 


E( error) 


=  Fn/3 


,-b 


(  A38) 


Even  if  the  worst  case  occurs  and  errors  do  not  cancel  at 
all  in  taking  FFT  and  IFFT,  then  the  total  convolution  error 


E ( error) 


JN/3  2  100% 


(  A39) 


To  give  a  numerical  example,  for  the  IBM  mainframe 


computers,  b=23.  For  N=100000, 


E ( error ) 


,100000/3  x  2~22  x  100X 


(  A40) 


=  .  004  X 

Therefore,  if  we  check  the  answer  by  comparing  V  with  V 

3 

computed  from  ( A 35 )  by  using  FFT  and  IFFT ,  the  numerical 

— * 

error  in  V  would  be  in  the  fifth  precision  position  i.e. 

3 

insigni f leant . 


Errors  due  to  ( b)  will  be  small,  provided  that  the 
problem  is  well  behaved  (i.e.  the  condition  number  of  the 
impedance  matrix  [Z]  is  small)  and  that  the  iterations  are 
carried  out  far  enough  so  that  the  error  in  V  is  small.  In 
a  practical  problem  like  the  antenna  array  problem,  the 
change  in  current  for  each  iteration  drops  off  very  sharply, 
indicating  that  the  actual  error  in  the  current  is  probably 
less  than  the  last  change.  The  following  qualitative 
argument  can  be  given  in  support  of  the  above  claim.  From 
(  A22)  , 

7n+1  =  [Q]  V  +  [R]  Jn  (A41  ) 

Jn  =  C Q ]  V  +  [R]  7n_1  (A42) 

Subtracting  (A42)  from  (A41),  we  get 

.  cm  !„  U43> 

?  *♦  ■*  th 

where  d  s  J  -J  ,  is  the  change  in  current  at  n 
n  n  n-l 

iteration*  Therefore 
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IV.  COMPUTER  PROGRAMS  AND  SUBROUTINES 

The  computer  programs  and  subroutines  given  in  this 
section  are  written  to  verify  that  DCM  works  properly  and 
to  measure  the  number  of  iterations  needed  for  some  sample 
problems.  No  attempts  have  been  made  to  optimize  the  computer 
code.  In  fact,  the  fast  fourier  transform  and  inverse  fast 
fourier  transform  (FFT  and  IFFT)  routines  given  are  for  2N 
points.  This  means  that  more  points  than  are  strictly 
necessary  has  to  be  taken.  However,  even  with  the  relatively 
unoptimized  code,  DCM  proves  to  be  faster  than  other  techniques 
for  large  problems. 

Subroutines  FFT  and  IFFT  are  for  one  dimensional  FFT 
and  IFFT.  Subroutines  TWODF  and  ITWODF  are  for  two  dimensional 
FFT  and  IFFT.  The  main  program  segment  starting  on  page  47, 
solves  the  one  dimensional  convolution  equation.  This  program 
is  written  to  be  able  to  solve  problems  with  gaps.  Subroutine 
SOLVE  solves  the  two  dimensional  convolution  equation.  It  is 
not  written  to  solve  the  problems  with  gaps,  however. 


SOBBOOTIME  FFT(Z,B,B) 

IBTEGEB  I,J,N,8,LE,LE1,HV2,H81,K 
COBPLEZ  1(4096), 0,8,1 
BEAL  El 

PI=3. 1415926535 
00  20  1*1,1! 

LE=2**(8*1-L) 

LE1=LE/2 

0*(1.0,0.0) 

8*CBPLZ  (COS  (Pl/f LOAT  (LEI))  ,-SIB  (PI/ FLOAT  (LE 1)  )  ) 
DO  20  J*1,LE1 
00  10  I=J,B,L£ 

IP=I*LB1 
T*Z(I)  +  Z(IP) 

Z(IP)  =  (Z(I)-Z(IP)  )*0 
I  (I)  =T 
10  COBTIBOE 
0*0*8 

20  COHTIMUE 
■f 2*8/2 
881=8-1 
J*1 

DO  40  1=1,881 
IF(I.GE.J)  GOTO  25 
T*Z  (J) 

Z  (J)  =  Z  (1) 

Z(I)=T 

25  K*BV2 

26  IF(K.GE.J)  GOTO  30 
J= J-K 

K*K/2 
GOTO  26 
30  J=J*K 
40  C0BTI8UE 
BET0B8 
EBD 

SOBBOOTI8E  IFFT(Z,8,8) 

I8TEGEB  I,J,8,B,LE,LE1,8T2,881,K 
COBPLEZ  1(4096) ,0,8,1 
Bill  PI 

PI* 3, 1415926535 
DO  20  1*1,8 
LE»2**(8*1-L) 

LI1*LE/2 

0*(1.0,0.0) 

B*CHPII  (COS  (PI/FLOAT  (LEI)  )  ,  SIB  (PI/FLOAT  (IB1)  )  ) 

DO  20  J*1,LEl 

DO  10  I*J,8,  LB 

IP*I+LE1 

T*I  (I)  +1  (IP) 

Z(IP)*(Z(I)-Z(IP))*0 
Z  (I)*T 
10  C08TI80E 
0*0*8 

20  C08TI80E 
Bf2*B/2 
881*8-1 
J*1 

DO  40  1*1,881 
IF(I.GE.J)  GOTO  25 
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T*X  (J) 

X  (J)  =  X  (I) 

X|I)=1 

25  K=MV2 

26  IF(K.GE.J)  GOTO  30 
J*J— K 

K-K/2 
GOTO  26 
30  J=J*K 
40  CONTINUE 

00  50  1=1,11 
X  (I)  =  X  (I) /FLOAT  |M) 

50  CONTINUE 
BETUBN 
END 

C**  DISCHETE  CONVOLUTION  METHOD  10  SOLVE  BOTH  TOEPLITZ  END 

C**  NON  TOEPLITZ  NATBIX  EQUATIONS  B¥  (N  LCG  N  PBOCESS)  ITEBATION 

IHTBGEB  N,  N 

COMPLEX  X ( 1 366) ,T(4096) ,CZEBC, CUB  (4096) ,1(4096) ,TSTORE 
INTEGEB  STABT  (11),FINISH(10) ,ST,FI,BEVS 
C2EBO=(0. 0,0.0) 

1  DO  10  1=1 , 1024 
A  (I) =CZEBO 
CUB (I) =CZ  EHO 
HI)  =CZEBO 
10  CONTINUE 

BEAD(1, 100) N,H, NS,IFLAG,IBEG,BEVS 
BBITE  (3,200) N,H, MS 
IF(IPLAG.NE.O)  GO  TO  11 
BEAD  (  1,300)  (X  (I), 1  =  1, NS) 

MBITS  (3,  400)  (X  (I)  ,1  =  1,  IS) 

BEAD  (1,300)  (T  (I)  ,1=  1 ,  NS) 

HBITE  (3,  600)  (T  (I)  ,1  =  1, NS) 

GO  TO  13 

1 1  CONTINUE 
NSBT2=MS/2 

OPEN ( 0NIT=2 1 ,FILE=»  ABDATA. DAT ' ) 

BEAD  (21,101)  (T  (I), 1=1, NS) 

BEAD  (21,101)  (X  (I)  ,1=1,  NS) 

101  FOBHAT (5E 1 4*  7) 

IF(BEVS.EQ.O)  GO  TO  13 
DO  12  I=1,NSBT2 
TSTOBE-T (I) 

IPTB=MS-I«-1 
T  (I)  *T  (IPTB) 

T  (IPTfi) =TSTOBE 

12  CONTINUE 

13  CONTINUE 

IF (IBBG.EQ.O)  GO  TO  14 

BEAD f  1, 100)  (STABT  (I)  .FINISH(I)  ,I*1,IBEG) 

IN  CONTINUE 
NS2*NS*MS 
NSN1=NS-  1 
DO  15  1*1, NS 
J*X+NSfl 1 
A  (J)*X  (I) 

T  (NS2-I)  =T  (I) 

15  CONTINUE 

CALL  FFT (I,N,B) 

I COO NT- 1 
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IF  (IHEG.EQ. 0)  GO  TG  18 
00  16  1=1 , IBEG 
ST-STABT (I) 

FI-FINISH (I) 

DO  16  J-ST ,FI 
A(J)=CZEBO 
16  CONTINUE 
18  CONTINUE 
20  CONTINUE 

CALL  FFT(A,N,N) 

00  30  1*1,  B 
COB (I) =A (I) /T  (I) 

30  CONTINUE 

CALL  IFFT  (CUR, N,H) 

■BITE  (3,1000) ICCUNT 
DO  40  I-NS*1,N 
COB (I) =CZEBO 
40  CONTINUE 

IFIIBEG.EO.O)  GO  TO  46 
DO  45  1*1, IB EG 
ST-STABT (I) 

FI-FINISH (I) 

DO  45  J-ST, FI 
COB  (J) =CZEBO 

45  CONTINUE 

46  CONTINUE 

■BITE  (3,600)  (CUB  (I)  ,1*1, NS) 

CALL  FFT (CUB, N,E) 

DO  50  1*1,  ■ 

A  (I)  -CUB  (I)  *T  (I) 

50  CONTINUE 

CALL  IFFT  (A, N,H) 

■BITE  (3,500)  (All)  ,I-BS,BS»NSHl) 

COHEBB-0.  0 

EBAX*0.0 

IF (IBEG. BE. 0)  GO  TO  52 
DO  51  1*1, NS 
J«I*NSH1 
TSTOBE-X  (I) 

EBBOB*CABS  (A  (J) -TSTOBE)  /CABS  (1ST0BE) 
A (J)*TSTOBB 

IF (EBBOfi.GT. EHAX)  E0AX-EBBOB 
CO0EBB-CU0EBB+EBBOE 

51  CONTINUE 
GO  TO  54 

52  CONTINUE 
ST*1 

DO  53  1*1, IBEG* 1 

IF  (I.  NE.  1 )  ST-FINISH(I-I)  *1 

FI*STABT(I)-1 

IF  (I.  EC-IBEG*1)FI*NS 

DO  53  J*ST,FI 

K-J+NS01 

TSTOBE*! (J) 

BBBOB*C AES  (A  (R)  -TSTOBE)  /CABS  (TSTOBE) 
IF(EBBOB.GT.EBAX)  E0AX-EBBOB 
C0BEBB»C0HBBB*EBBG8 
A (R) -TSTOBE 

53  CONTINUE 

54  COBT2NOE 
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EHAX=EH1X*  100-0 

COMER 8=100. 0*COHEBB/FLOAT (NS) 

■BITE  (3,900)  EHAX.CDMEBB 
MBITS  13.700) 

READ! 1. 100) IFLAG 
I?  (If  116.  BQ.  0)  60  TO  55 
C1L1  IFF!  (COB.M.fl) 

■HITE  (3.600)  (COB  (I)  ,1=1.  NS) 

■BITE  (3,500)  (1  (I).I=MS,MS»MSB1) 

IF  (1FL1G.  EQ.  1)  STOP 
■BITE (3, 1 100) 

60  TO  1 
55  COMIIMOE 

IC0UMT=IC00MT*1 
60  TO  20 

100  FOBHAT  (1010) 

200  FOBHAT  (IB  N  M  MS*. 317) 

300  FOBHAT  (10B0.0) 

400  FOBHAT  ( 1H  , *  EXCITATION  VECTOB*/(lH  .10E11.4)) 

500  FOBHAT (1H  .'RESULTANT  FIELD'/(1H  .10B11.4)) 

600  FOBHAT  (  1H  .  'CUBBENTS  '/ ( 1H  .10E11-4)) 

700  FOBHAT  <1H  ,»COMT INUE  ITERATIONS?  0  FOR  YES,', 

$  •  1  FOB  MO  AMO  BETOBM*) 

800  FOBHAT  (  1H  ,  *  GBEEMMS  f  ONCTIOH*  /  (1 OE 1 1.  4)  ) 

900  FOBHAT  (IB  ,'HAX  FIELD  EHHOH= • , E 10. 3 , •* */ 
t  IB  , 'AVERAGE  EBBOB  =  • , E 10-3 , •*  • ) 

1000  FOBHAT  (1H  ,//1H  , ‘ITEBATIOM  MOHBEB • ,17//) 

1100  FOBHAT  (IB  ,//1H  , ' **********MBXI  PBCBIFH**********//) 
EMD 


SOBBOOTIME  PFT  (X, fl, STABT, STEP) 
COMPLEX  X  (16384) ,0,1,1 
IMTEGEB  ST ART, STEP, SPIFF 
M*2**B 

SDIFF=ST1P-STABT 
HV2  =N/2*STEP 
Mfll  =  (N-2)*STBP*STABT 
M  =  (N-1)*STBP+ STABT 
J  =STABT 
DO  8  I=ST ABT ,MH 1 , STEP 
IF(I.GE.J)  GO  TO  5 
T  =  X(J) 

X|J)  =  X  (I) 

X(I)  =T 

5  K  =  MV2 

6  IF(K-SDIFF.GE.J)  GO  TO  7 
J  =  J-K 

K  =K/2 
GO  TO  6 

7  J  =J*K 

8  COMTINOE 

PI  =3.14159265358979 

DO  20  L=  1 , 8 

LB  =2**L 

LE1=LE/2 

LSTEP=LE 1 *STEP 

U  =(1.0, 0.0) 

AMGLE=PI/FLOAT  (LEI) 

H  =CMPLX  (COS  (AMGLE)  ,  -SIM  (ANGLE)  ) 

LEI  =LSTEP*STABf -STEP 
LE  =LE*STEP 

DO  20  J=STABI, LEI, STEP 
DO  10  I=J, N,LE 
IP  =I*LSTEP 

T  =X(IP)*0 

X(IP)  =  X  (I)  -T 
X  (I)  =X(I)*T 

10  COMTIMUE 
0=0*11 

20  COMTINOE 
BETOBM 
BMD 

SOBBOOTIME  THODF  (X,M,B,L,BR,LB) 
COMPLEX  X  (16384) 

IMTEGEB  STABT, STEP 

SIABT=1 

STEP  =1 

DO  10  1=1, fl 

CALL  FFT(X, LB, STABT, STEP) 
STABT=STABT+L 
10  COMTINOE 
STEP  *L 
DO  20  1*1, L 
STABT* I 

CALL  FFT(X,flfl, STABT, STEP) 

20  COMTINOE 
BETOBM 
EID 

SOBBOOTIME  IFFT  (X,fl, STABT, STEP) 
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COflPLEX  X  (16384), 0,B,T 
IBTEGER  START, STEP,SDIFF 

SDIFF*STEP— START 
Mf 2  *B/2«STBP 
ESI  *(B-2)*STEP*STABT 
BEXP  * <B- 1) *STEP+SIABT 
J  -START 
DO  8  I*START,IH1, STEP 
IF  |I.  GE.  J)  GO  TO  S 
T  *X(J) 

I(J)  *X|I) 

HI)  *1 

5  K  =H?2 

6  IF  (K-SDIPF.GE. J)  GO  TO  7 
J  =J-I 

K  =K/2 
GO  TO  6 

7  <3  *J+K 

8  COBTIBOE 

PI  =3.14159265358979 
DO  20  1=1,  fl 
LB  =2**L 

LEl*LB/2 
LSTEP=LE 1*STBP 
0  =11.0,0.0) 

ABGLE=PI/FLOAT|LE1) 

B  =CflPLX  (COS  IABGLE)  ,  SIB  ( ABGLE)  ) 

LEI  =LSTEP* START-STEP 
LB  *LI*STEP 

DO  20  J=START, LEI , STEP 
DO  10  I=J, BEXP, LB 
IP  =I*LSTEP 

T  =X(IP)*0 

X  (IP)  *X  (I)  -T 
X(I)  =X  (I)  +T 

10  COBTIBUE 
0*0*  B 

20  COBTIBOE 
BETOBB 
BBD 

SOBBOOTIBE  IT80DF  (X,B,fl,L,Bfl,LB) 
COflPLEX  X(  16384) 

IBTEGEB  START, STEP 
START* 1 
STEP  *1 
DO  10  1*1, fl 

CALL  IPFT  (X,Lfl, START, STEP) 
SIABI»START*L 
10  COBTIBOE 
STEP  *L 
DO  20  1*1, L 
SIABT*I 

CALL  IFFT (X,RH,STABT,S1BF) 

20  COBTIBOE 
PB*FLOAT(B) 

DO  30  I*  1, B 
X(I)*X(I)/FB 
30  COBTIBOE 
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SOBBOUTINB  SOLVE(A,fi,LO,BO,BO,LB,flR,L,B,N,FXCITE) 

BOOTIME  TO  SOLTI  TBE  BATBIX  EQUATION  A  X  =  B 
LOGICAL  PICITE 

COMPLEX  CTEHP,CZBBC,  A  (  16384)  ,1116384)  ,V (16384)  ,B(1849)  ,1(1849) 
IBTB6BB  FLAG  1 , FLAG 2, CO UNI 
BBAL  CHNAYG,CHNNAX, CHANGE 
CSBBO-  (0.0, 0.0) 

ZIBOIZE  I  AMO  V  (EXPENDED  E) 

DO  10  1=1, N 
T  (I)*CZ£BO 
COVTINOE 
00  20  Is  1 ,  MO 
Y  (I) =CZEBO 
COBTIBOE 

IP  (.HOT. PXCITE)  GO  TO  65 
00  40  I»  1,(10 
IPTH= (I- 1 ) *L 
JPTB=IPTB*LO*LO- 1 
DO  30  J=IPTB*1,IPTB*LO-1 
A  (JPTB)*A  (J) 

JPTB=JPTB-1 

CONTINUE 

CONTINUE 

FILL  UP  A  ABBAY  AND  V  ABBAY 
DO  60  1=1, HO- 1 
IPTH*  (1-1 )  *L 
JPTB* (BO*BO-I- 1) *L 
DO  50  J= 1 ,LO*LO- 1 
A  (JPTB+J)  =A  (IPTB* J) 

CONTI IUE 
CONTINUE 

CALL  THOOF (A,M,fl,L, Nfl,Lfl) 

COUNT*0 
JPTB=  1 

COONT*COUMI+1 
DO  90  1*1,110 
IPTB»L*(NO-2*I)*IC 
DO  80  J*  1 ,  LO 
T (IPTB) *B  (JPTB) 

JPTB* JPTB* 1 
IPTB*IPTB* 1 
CCBTINOE 
CONTINUE 

FIND  T  TEAMS FOBBED  AND  COMPUTE  X  TBANSFOBHBD 
CALL  IHODF  (¥,M,H,L,Bfl,LH) 

DO  100  I* 1, B 
X  (I)  *T  (I)  /A  (I) 

CONTINUE 

GET  X  FBOB  X  TfiANSFOBHED 
CALL  ITBODF(X,M,B,l,RN,LB) 

TBONCATE  X  AND  SATE  X  AFTBB  COBPUTIMG  TBB  CONTBBGENCB  CBIT1BION 

CHNATG*0. 0 

CBNNAX*0. 0 

DO  120  I*  1,N 

IPT1-1/L 

JPTB*I-IPTB*L 

If (JPTB. LE.LO.AND. JPTB. BE.O.AND.IPTB.IT.BC)  GO  TO  110 

XID-CZEBC 

60  TO  120 

IPTB*IPTB*LO*JPTB 
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CTEMPSX |I) 

CHANG  B=CABS (CTEMP-Y (IPTB) )/CABS (CTEMP) 

I(IPTH)=CTEMP 
CBMAVG=CHMAVG*CHANGE 
IP  (CHANGE. GT.CHNMAX)  CHNHAX-CBANGE 
120  COMTINUB 

CHNAVG3  (CHNAVG*  100.0) /FLOAT (NO) 

CBNHAX=CHNHAX*1G0.0 
MBITS (3,1200)  CHNAVG, CHNMAX, COUNT 
C  FIND  THB  TBANSFOBH  OF  TBOBCATEO  X 

CALL  THODF  (X,N,M,1,MH,LM) 

C  COMPOTE  V  TBAMSFOBHED 

00  130  1=1, N 
VfI)=A(I)*X(I) 

130  CONTINUE 

C  GET  V  FBOfl  V  TBAMSFOBHED 

CALL  ITHODF (V, N,fl,l, N9, LB) 

C  COMPUTE  THB  BBBOB  CBITBBION 

CHNAVG=0. 0 
CHNHAX=0.  0 
JPTB= 1 

DO  150  1=1, NO 
IPTB=I* (MO-2*I) *10 
DO  140  J= 1,LO 
CTENP=B  (JPTB) 

CHANGE=CABS  (CTEMP- V  (IPTB)  )  /CABS  (CTBHP) 

CHNAVG=CHNAVG*CHANG£ 

IF (CHANGE. GT.CHBMAX)  CHNMAX=CHANGB 
IPTH=IPTH* 1 
J PTH = JPTB* 1 
140  CONTINUE 
150  CONTINUE 

C  ASK  HHETHEB  OB  NOT  TO  STOP  AFTEB  BEPOBTING  X  FIELD  EBBOB 

CHNAVG=  (CHNAVG*  100.0)  /FLOAT  (NO) 

CHNHAX=CHNMAX* 100.0 
MBITS (3, 1300)CHN BAX, CHNAVG 
BEAD  (1,1100) FLAGl 
IF(FLAGI.EQ.O)  GO  TO  70 
MBITE  (3,  1400)  (I  (I)  ,  1=  1,NO) 

C  ASK  IF  FIELD  SHOULD  BE  PBINTED  CUT  ALSO 

MBITS  (3, 1500) 

BEAD  (1,1 100)  FLAG2 
IF  (FLAG2.  NE.  0)  BETUBM 
MBITE  (3,  1600)  (V  (I)  ,  Is  1 ,  N) 

BETUBM 

1000  FOBMAT ( 10E0. 0) 

1100  FOBMAT  (410) 

1200  FOBMAT (Ifl  ,'A¥G  CUBBENT  CHANGE3 ' ,E 1 4. 7, '  X'/ 

$  1H  , 'MAX  CUBBENT  CHANGE3* ,E 14. 7, •  X*/ 

$  IB  , 'AFTEB*, 14,*  I1EBATICNS*/) 

1300  FOBMAT  (1H  ,*MAX  FIELD  EBBOB  3  *,E15.7,»  X*/ 

$  1H  , * A VG  FIBLD  EBBOB  3  *,E15.7,*  X*/ 

$  IB  , 'CONTINUE  IIEB ATIONS?  0  FOB  IBS,  1  FOB  NO,  AND  BETUBN*/) 

1400  FOBMAT  (1H  , 'CUBBSNTS'// (IB  ,10X11.4)) 

1500  FOBMAT  (IB  ,*PBINT  FIELDS?  0  FCB  1ES,  1  FOB  NO,  TBII  BETUBN'/) 
1600  FOBMAT  (1H  ,'BESULTANX  FIELDS'//(1B  ,10X11.4)) 

END 
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